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Introduction 


In this paper we present a new form of the collocation method that 
allows one to find very accurate solutions to time marching 
problems without the unwelcome appearance of Gibb’s 
phenomenon oscillations. The basic method is applicable to any 
partial differential equation whose solution is a continuous, albeit 
possibly rapidly vaiying function. Discontinuous functions are 
dealt with by replacing the function in a small neighborhood of the 
discontinuity with a spline that smoothly connects the function 
segments on either side of the discontinuity. This will be 
demonstrated when the solution to the inviscid Burgers equation is 
discussed. 

As is well known it is possible to represent a smooth function with 
great precision if the representation is made in a basis set of 
solutions of a singular Sturm-Liouville problem. However if the 
function varies rapidly unreal oscillations in the representation may 
occur unless a very large number of terms are used. We consider 
an alternative approach of dividing the domain of the function into 
an arbitrary number of subdomains. In each subdomain the 
function segment will be simple enough to be represented in a low 
order expansion. Differential equations will be solved separately 
in each subdomain. Truncation errors inevitably follow from using 
a finite order expansion and these errors manifest themselves most 
noticeably as discontinuities across subdomain boundaries in high 
order derivatives. The tau method will be used to establish 
continuity of the highest order derivative that occurs in the partial 
differential equation under consideration. As the function evolves 
in time the subdomain boundaries will be permitted to move. In 
this way the growth of instabilities and Gibb’s phenomenon 
ringing will be avoided. 



Subdomains 


The ideal set of functions to span a finite subregion of a possibly 
very large domain is the set of Chebyshev polynomials. If a 
function segment is prescribed in the subdomain x^ < x < x^ of 

the overall domain x,<x<x and the mapping x'= a k x + b k such 
that x'= -1 when x = x mm and x'= +1 whenx = x^ is performed, 
then the function at x may be represented by a Chebyshev 
polynomial expansion at x’ : 

f<x) = "tXT,(x') (1) 

l~0 

The index k identifies the particular subdomain in which x is 
found. We evaluate the function at the collocation points: 

( 2 ) 

where /, =-i + N , and 1 <i<N, with this numbering scheme 
x\ =-l and x N =1. Then 

/(^) = 'f UX(x,) (3) 

a k 

or in matrix form: 

YJA ‘=f,‘ (4) 

if N is small the expansion coefficients may be found quickly and 
accurately by matrix techniques (either inversion or LU 
decomposition) 



If the function segment is well represented by N polynomials then 
the highest order polynomials should be associated with very small 
coefficients. Using a simple algorithm discussed below a set of 
contiguous subdomains can always be found in which the last four 
expansion coefficients are very small: 


A k N _ 3 <£, A k N _ 2 <£ , A*_, <£, A k N <£ , £«\ 


In Figure 1 is a superposition of the function y = e~* and a 
Chebyshev representation made in three subdomains with eighteen 
polynomials in each subdomain: 



Figure 1 : A superposition of a Gaussian and its Chebyshev 
representation. The curves are virtually indistinguishable 

Figures (2)-(4) show each subdomain separately. 



Figure 2. The Chebyshev representation in the first subdomain 



Figure 3. The Chebyshev representation in the second subdomain 
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Figure 4. The Chebyshev representation in the third subdomain 



Figures (5) and (6) demonstrate the considerable accuracy possible 
with this method. Figure (5) is a superposition of the nearly 
singular function: 


-2sinh(x) 
cosh(x) - e" 0001 


( 5 ) 


and its Chebyshev representation. Figure (6) is a superposition of 
the second derivative of the function and the second derivative of 
its Chebyshev representation: 




Figure 5. A nearly singular function and its Chebyshev 

representation 



Figure 6. A superposition of the second derivative of the 
function and the second derivative of its Chebyshev 
representation The representation of the function and of its 
rapidly varying second derivative is excellent and there is no 
sign of Gibb’s oscillations. 



To form a strategy for finding subdomain boundaries we consider 
the problem of solving the Burgers equation 

fyM _ ( gMfjO (6) 

a/ ^ ' dx dx 2 

y(x,t 0 ) is the function shown in figure (5). Again the function 
will be treated as a set of function segments in contiguous 
subdomains each represented by a low order Chebyshev 
polynomial expansion. We will use the tau method to find 
expansion coefficients in each subdomain. The tau method adds to 
the N order expansion m additional terms of order N+l to N+m, 
which impose m boundary conditions on the polynomial 
representation, and it will be used to ameliorate the effects of the 
most serious source of error, which is the truncation error that 
occurs in a low order polynomial expansion. The tau method will 
be used to enforce continuity of the highest order derivative that 
occurs in the partial differential equation across subdomain 
boundaries. 

Assume that the subdomain labeled k begins aix a and ends at x b . 

In this subdomain x is mapped into x so that - 1 < x < 1 , 
andx 5= a k x + b k . 

We consider a tau method expansion defined by: 

I 'A,T,{x)= f{x), \<i<N, (7) 

/=0 


and 

!)=/'&) 

af£V/(+l) = 

/=0 


(7a) 

(7b) 



( note that 


( dx') 


dx 2 dx 


dx 


-a 


J 


1 dx 


)• 


Eqs. (7) form a set of N+2 equations that produce a set of 
expansion coefficients which represent the function exactly at N 
collocation points and correctly reproduce the second derivative of 
the function at the endpoints of the subdomain. 

These equations may be combined into the single matrix equation: 


XS„A,=f>, 

S,=T,(x,), l<i<N , 0</ <iV + l 

s„ M =r(-i). 0 </<at+i 

Oil<N + l 

and 



1 < z < 


( 8 ) 

(8a) 

(8b) 

(8c) 


(8d) 

(8e) 

(8f) 


At the inner boundary of the innermost subdomain and at the outer 
boundaiy of the outermost subdomain second derivatives need not 
be specified because there are no neighboring subdomains to 
match boundary derivatives there. Accordingly for the first 
subdomain we consider a set of matrices similar to those in eqs. (8) 
except that the value of the second derivative is needed only at the 
outer edge: 


2X4=*, 


(9) 


where 


X,=T{x), 1 <i<N, 0 <1<N 

x nm =t;(+\), o<i<n 

and 


f 


h,=y 


x\ -b, 







\<i<N 


(9a) 

(9b) 


(9c) 

(9d) 


Solving eqs. (9) for the^, produces a set of N+l coefficients 
0 <1 <N, that exactly reproduce the function at N collocation 
points, and the correct second derivative at the outer boundary of 
the innermost subdomain. We include the definition: 


X., - 0 (9e) 

to create a set of N+2 coefficients as in eqs.(8). 


Similarly for the last subdomain: 

2X4=*, ( 10 > 


where 


Y„~T,(x), \<i<N, 0 <1<N, 


(10a) 


r«,-JT(- 1). o<i <n , 


(10b) 


and 



\ 

y 



= o 


l<z<W 


(10c) 

(lOd) 

(lOe) 


Finally there exists the possibility that one subdomain is 
sufficient in which case the subdomain is identical with 
the complete computational domain and no boundary 
matching needs to be done. Then 


YJA=h 


(ii) 


where 


f 


h, =y 


x'-b t 


V 


o h 


1 <i<N 


(11a) 


Subdomain Boimdaries 

If the domain of y(x) is x a < x <x b then letx = x a and x f = x b and 

test the possibility that one subdomain is sufficient. To do this 
simply map the domain into that appropriate for Chebyshev 
polynomials: x a — >• -l,x b -> +1, evaluate y(x) at the collocation 
points, and invert equation (11): 

4 = ZZ; 1 / for 1 < i < N, and 0 < / < N -l 

i 


( 12 ) 



to find the Chebyshev coefficients and we define 

4=o, 4„ =o, 


(12a) 


to create a set of N+2 coefficients. If the representation is a good 
one then the coefficients of the highest order polynomials will be 
very small. We apply a test to these high order coefficients: 


(1.) Let = 1 , then. 

Find the coefficient whose magnitude, \A. is largest. If this 
quantity is larger than then let A imx = A j . 


(2.) Compare the highest order coefficients with 4^. They must 
satisfy the smallness test: 




<€, 


N -4<m<N + 1 , £ « 1 , 


(13) 


If the smallness test is passed then one subdomain is enough. 
Otherwise with x = x and x = x . : 

i a jo 

(1 .) Let x = x t , and redefine Xf so that 

x f =c l x i +c 2 x f where c, +c 2 =1. (14) 

Then use eqs.(9) (This set of equations is used because x. is the 
innermost point in the computational domain), to find the N+2 



expansion coefficients and again apply the smallness test. If the 
test is failed then repeat (1 .) until it is passed. This is the 
innermost subdomain extending from x]' = x, to x ] ‘ = x f . Then: 

(2.) Letx. = x x ‘ and x f = x b . Use eqs. (10) to find the expansion 
coefficients (This set of equations is used because x f is the 

outermost point in the computational domain), and apply the 
smallness test. If it is passed then two subdomains are sufficient. 
Otherwise: 

(3.) Let x f - x b and 

x i =c 2 x i +c 1 x f (15) 

Use eqs. (10) to find the expansion coefficients and again apply the 
smallness test. If the test is failed then repeat step (3.) until a 
sufficiently small subdomain is found. This is the outermost 
subdomain extending from*' 0 = x, tox^° = x b . 

(4.) Return to the inner boundary. Letx, = x u f and x f = x)° . Use 

eqs. (8) to find the expansion coefficients and apply the smallness 
test. If the test is passed then three subdomains are sufficient. If 
not then: 

(5.) Let x. = x, andx f = c,x. + c 2 x f . Use eqs. (8) to find the 

expansion coefficients, and apply the smallness test. If the 
expansion coefficients fail this test then repeat step (5.) until a 
sufficiently small subdomain is found. This is the second 
innermost subdomain. It extends from x, 2 ' = x" to x 2i = x f . 

(6.) Return to the outer boundaiy. Letx, = x 2 'andx / = x)° Use eqs. 

(8) to find the expansion coefficients and apply the smallness test. 
If the expansion coefficients pass this test then four subdomains 
are sufficient. If not then: 

(7.) Let x f = x'° and x, = c 2 x, + c l x / . Use eqs. (8) to find the 

expansion coefficients and again apply the smallness test. If it is 
failed then repeat step (7.) until a sufficiently small subdomain is 



found. This is the second outermost subdomain extending 
from x ; 2 ° = x. tox 2 ° - x ]° . 

(8.) Repeat this search pattern, (i.e. next find* 3 ' andx 3 ' , and 
then x 3 ° and x 3 °, etc.), until a complete set of subdomains is found. 


Partial Differential Equations 

The solution of the partial differential equation in each subdomain 
is found through integration over a small time interval followed by 
a corrector step to match derivatives at subdomain boundaries. 


Returning to the Burgers equation as illustration, we see that the 
right hand side of eq. (6) is: 

( 16 ) 

y:(t) = a l 'fT,A!(t) (17) 

y:‘(‘)=A'T-A;(t) (18) 

/=0 

in each subdomain. 

Ignoring boundary conditions for the moment; the differential 
equation in each subdomain becomes 



dAjjt) 

dt 




(19) 


Boundary equations are established by specifying g’, g ! °* and by 
averaging the P.D.E. operator at the outer boundary of one 
subdomain and the inner boundary of the neighboring subdomain. 



( 20 ) 


Q m - (t)y? i‘)+yT' (Q-y! (<)y($ + y? (0 

*rw-*r« (2D 

then 

~~=v g :(t) (22) 

is an ordinary differential equation for N expansion coefficients 
which may be solved using standard methods such as Runge-Kutta 
integration. The solutions to eq. (22) will maintain continuity of 
the function across subdomain boundaries but the derivatives are 
unconstrained and can become noticeably discontinuous. To 
impose continuity of the highest order derivative that occurs in the 
P.D.E. at subdomain boundaries we invoke the following 
corrector: 


We begin with eqn (22) which provides an initial expression for 


dA,‘(r) 

dt 

point: 


, and this is used to calculate 


< fr,‘(0 

dt 


at each collocation 


dy k ,(t) dA‘(t) 
dt ^ * dt 


(23) 


and the change in time of the highest order derivative at 
each boundary point: 



dt 


dt 


( 24 ) 



(25) 


dy;-‘ _ df(t) 

®kZu*\l 


dt 


dt 


J ' «,£ J ' n,k-l 

and — are the same quantities evaluated at the common 
dt dt 

boundary of neighboring subdomains. To enforce equality the 
corrector step replaces these quantities by their averages: 


dy"‘ , dy" K 


n 9 k-\ 


dy" k dt 
dt 


+ 


dt 


( 26 ) 


dyT _ dy 


n,k 


dt 


dt 


( 27 ) 


dy r 1 _ dy- 1 


dt 


dt 


( 28 ) 


and uses these averages in matrix operations similar to those 
defined by eqs. (8). For example in an interior subdomain we have 


i dt 


( 29 ) 


where 


h k = 


<W(t) 

dt 


1 <i<N 


h* = 


N < fl 


dy;\t) 1 
dt a ; 


,, «W(0 1 

dt a : 


( 30 ) 

( 31 ) 


(32) 



Then 


^^ = 1 IS-'h'it), 0<1<N + 1 (33) 

dt « * 

is an ordinary differential equation for the rate of change of the 
expansion coefficients which maintains continuity of the function 
and its highest order derivative at the common boundary of 
subdomain k and subdomain k+1, and at the common boundary of 
subdomain k and k-1 . 


Moving Boundaries 


Repeated integration of equation (33) results in the growth of 
discontinuities in high order derivatives at subdomain boundaries. 
To stop the growth of these discontinuities it is necessary to halt 
the process from time to time and find a new set of subdomains 
whose boundaries are different from those of the present set. 

There is no unique rule for determining when this should be done, 
but one effective rule follows from monitoring the growth of the 
high order polynomial coefficients. If originally the function in 
each subdomain is described by the polynomial coefficients A° r,gma ‘ 


the high order terms are defined to be those for which 



<8 


for all / > L . As time goes on we test for the condition 



^ original 


>8 


for / > L 


(34) 


When this condition is met the high order terms are growing and 
new boundaries must be found. To do this, simply stop the 
integration routine, treat the present state as the initial condition 
and reapply the algorithm for finding subdomains discussed above 



except that a search from Xinitiai to x^i should alternate with a 
search from x fina i to x^t^i . 

In summary solving P.D.E.’s using this localized collocation 
method is a three-step process: 

(1 .) Find a set of subdomains in which a function is represented 
by a small set of Chebyshev polynomials. 

(2) Solve the O.D.E. eqn (33) in each subdomain. 

(3) Test for the growth of high order coefficients with eqn (34). 

If eqn (34) is satisfied then reinitialize the problem by returning to 
step (1). 


Examples 

The solution to Burgers equation without boundary conditions and 
with 


y(x,t = .0001) = 2sinh(x) as initial state is: 

cosh(x) - e 

y(x,t)= - 2sinh(Y) - (35) 

V cosh(x) - e 


A superposition of numerical and analytic solutions at various 
times is plotted in figure (7). The numerical solution has the 
boundaiy conditions that y{—\ y t ) and y(l,t) are constant in time 
The numerical solution is indistinguishable from the analytic 
except after a long time when the difference in boundaiy 
conditions becomes important. This is shown in figure 8. 




Figure 7 Solution to Burgers Equation. Each curve is 
a superposition of numerical and analytic solutions. 



Figure 8 Boundary Effects. The thick line is the analytic 
solution with no boundary conditions, the thin line is the 
numerical solution with no change allowed at the 

boundaries. 



A more interesting example is the Kortweg -Devries equation. 


dt dc 


(l+yfcO) 


dc* 


As the initial state we use: 
(=o ) = - 1 2 cosh(x) ” 2 


This equation has the soliton solution 
y(x, t) = -1 2 cosh(x - 50 " 2 
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(36) 


(37) 


(38) 


40.00 


Figure 9 Each curve is a superposition of the numerical and 
analytic solutions of the Kortweg-deVries equation at different 

moments in time. 


The wave equation is second order in time: 


d 2 ,y(x,0 _ d 2 y(x,t) 
dt 2 ~ dx 2 


( 39 ) 


with the initial state 


y(*,t t __ 0 ) = cos 


/ V“ 

' 7DC ) 




( 40 ) 


The figure below shows the solution at different moments 
in time. 



Figure 10. Initial pulse dividing in half and traveling in 

opposite directions 



Figure 1 1 Reflected wave which reassembles at the center 


Finally, to demonstrate the effect of the occurrence of a 
singularity consider the inviscid Burgers equation 


= (42) 

dt dx 

with 

y(x , t 0 ) = cos(^x) 4 , - 1 < x < 1 , (43) 

as the initial state. As the wave becomes steeper the size 
of the subdomains near the leading edge become smaller. 
When one falls below an arbitrary threshold in size a 
singular region is assumed to have been found, and the 
function in a small neighborhood of the singularity is 
replaced by a spline. The spline is a sixth order 
polynomial whose coefficients are determined by the 
value of the function and its first and second derivatives at 
the boundaries of the singular region. 




